pure subroutine profil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
!
! calling sequence:
! call profil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
!
! input parameters:
! iopt : integer flag, specifying whether the profile f(y) (iopt=0) or the profile g(x) (iopt=1)
! must be determined.
! tx : real array, length nx, which contains the position of the knots in the x-direction.
! nx : integer, giving the total number of knots in the x-direction
! ty : real array, length ny, which contains the position of the knots in the y-direction.
! ny : integer, giving the total number of knots in the y-direction
! c : real array, length (nx-kx-1)*(ny-ky-1), which contains the b-spline coefficients.
! kx,ky : integer values, giving the degrees of the spline.
! u : real value, specifying the requested profile.
! tx(kx+1)<=u<=tx(nx-kx), if iopt=0.
! ty(ky+1)<=u<=ty(ny-ky), if iopt=1.
! nu : on entry nu must specify the dimension of the array cu.
! nu >= ny if iopt=0, nu >= nx if iopt=1.
!
! output parameters:
! cu : real array of dimension (nu). on successful exit this array contains the b-spline
! ier : integer error flag
!
! restrictions:
! if iopt=0 : tx(kx+1) <= u <= tx(nx-kx), nu >=ny.
! if iopt=1 : ty(ky+1) <= u <= ty(ny-ky), nu >=nx.
!
! other subroutines required:
! fpbspl
!
! author :
! p.dierckx
! dept. computer science, k.u.leuven
! celestijnenlaan 200a, b-3001 heverlee, belgium.
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
! ..scalar arguments..
integer, intent(in) :: iopt,nx,ny,kx,ky,nu
integer, intent(out) :: ier
real(RKIND), intent(in) :: u
! ..array arguments..
real(RKIND), intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1))
real(RKIND), intent(out) :: cu(nu)
! ..local scalars..
integer :: i,kx1,ky1,l,l1,m0,nkx1,nky1
! ..local array
real(RKIND) :: h(MAX_ORDER+1)
! ..
! before starting computations a data check is made. if the input data
! are invalid control is immediately repassed to the calling program.
kx1 = kx+1
ky1 = ky+1
nkx1 = nx-kx1
nky1 = ny-ky1
ier = FITPACK_INPUT_ERROR
select case (iopt)
case (0)
if (nu<ny) return
if (u<tx(kx1) .or. u>tx(nkx1+1)) return
! the b-splinecoefficients of f(y) = s(u,y).
ier = FITPACK_OK
l = kx1
l1 = l+1
do while (u>=tx(l1) .and. l/=nkx1)
l = l1
l1 = l+1
end do
h = fpbspl(tx,nx,kx,u,l)
m0 = (l-kx1)*nky1+1
do i=1,nky1
cu(i) = dot_product(h(1:kx1),c(m0:m0+nky1*kx:nky1))
m0 = m0+1
end do
case (1)
if (nu<nx) return
if (u<ty(ky1) .or. u>ty(nky1+1)) return
! the b-splinecoefficients of g(x) = s(x,u).
ier = FITPACK_OK
l = ky1
l1 = l+1
do while (u>=ty(l1) .and. l/=nky1)
l = l1
l1 = l+1
end do
h = fpbspl(ty,ny,ky,u,l)
m0 = l-ky
do i=1,nkx1
cu(i) = dot_product(h(1:ky1),c(m0:m0+ky))
m0 = m0+nky1
end do
end select
return
end subroutine profil